Zebra

Package installation

Note: my local evalcast and covidcast installations may be slightly ahead of the remote versions due to outstanding PRs.

# remotes::install_github("cmu-delphi/covidcast", ref="main",
#                         subdir = "R-packages/covidcast")
# remotes::install_github("cmu-delphi/covidcast", ref = "evalcast-killcards",
#                          subdir = "R-packages/evalcast")
# remotes::install_github(repo="ryantibs/quantgen", subdir="R-package/quantgen")
# remotes::install_github("cmu-delphi/covidcast", ref = "main",
#                         subdir = "R-packages/modeltools")

Produce forecasts

Current model is quantile regression with all the same lags/covariates as zyzzyva.

• Signals:
  • “jhu-csse-confirmed_incidence_num”: lags 0, 1, 2, 3, 6, 9, …, 21.
  • “fb-survey-smoothed_hh_cmnty_cli”: lags 3, 10, 17, 24.
  • “doctor-visits-smoothed_cli”: lags 3, 10, 17, 24.
• We always have the forecast date as the second day of the epiweek (Monday).
• So we forecast 5, 12, 19, and 26 days ahead.
• No constraints on the predicted quantiles. But we sort after.
• No regularization. No CV.
• 1 model for the whole country.
• We use 28 days of training data.
• On each forecast date, we look at the most recent 83 days (as.numeric(forecast_dates[1] - start_day_zebra(forecast_dates[1]))) and keep the top 200 counties by total cases over the period.
• We apply corrections as in production.

# Quantile autoregression with 3 lags, or QAR3
 
corrector <- zookeeper::make_zyzzyva_corrector(
  zookeeper::default_county_params(data_source = "jhu-csse"))
zeb <- get_predictions(
  forecaster = quantgen_forecaster, 
  name_of_forecaster = "zebra",
  signals = tibble::tibble(
                      data_source = data_source, 
                      signal = data_signals,
                      start_day = list(start_day_zebra)),
  forecast_dates = forecast_dates, 
  as_of_override = function(x) lubridate::ymd(Sys.Date()),
  incidence_period = incidence_period, 
  ahead = ahead, 
  geo_type = "county", 
  signal_aggregation = "list", 
  geo_values = "*",
  apply_corrections = corrector,
  featurize = dav_maker_zebra,
  n = n, 
  tau = c(0.025, 0.100, 0.250, 0.500, 0.750, 0.900, 0.975),
  lags = lags, 
  lambda = 0, # Just do quantile regression 
  sort = sort, 
  nonneg = nonneg)

library(tidyr)
library(purrr)
library(ggplot2)
theme_set(theme_bw())
competition <- c("COVIDhub-ensemble","COVIDhub-baseline",
                 "CMU-TimeSeries", "Karlen-pypm")
submitted <- lapply(competition[1:3], get_covidhub_predictions, 
                    forecast_dates = forecast_dates, 
                    signal = "confirmed_incidence_num")
# submitted[[4]] <- get_covidhub_predictions("Karlen-pypm", 
#                                            forecast_dates = forecast_dates - 1,
#                                            signal = "confirmed_incidence_num") %>%
#  mutate(forecast_date = forecast_date + 1)
submitted <- bind_rows(submitted) %>% filter(ahead < 5)


# Some fixes to make comparable
qar_dc <-zeb %>% 
  mutate(quantile = as.numeric(quantile),
         incidence_period = "epiweek",
         ahead = ahead %/% 7 + 1) 

results <- evaluate_predictions(bind_rows(qar_dc, submitted),
                                backfill_buffer = 0,
                                geo_type = "county") %>%
  intersect_averagers(c("forecaster"), c("forecast_date", "geo_value"))

Overall AE, WIS, Coverage 80

NOTE: Results are based on the following numbers of common locations

results %>% group_by(forecast_date) %>% summarise(n_distinct(geo_value))

## # A tibble: 25 x 2
##    forecast_date `n_distinct(geo_value)`
##  * <date>                          <int>
##  1 2020-08-31                        164
##  2 2020-09-07                        165
##  3 2020-09-14                        166
##  4 2020-09-21                        167
##  5 2020-09-28                        168
##  6 2020-10-05                        167
##  7 2020-10-12                        166
##  8 2020-10-19                        166
##  9 2020-10-26                        168
## 10 2020-11-02                        168
## # … with 15 more rows

We compare the new forecaster to:

• COVIDhub-baseline
• COVIDhub-ensemble
• Our submission

Top line conclusions:

1. Overall performance approaches the ensemble.
2. By WIS, we are the top model for 1 weeks ahead, still good at 2-4 weeks ahead
3. Coverage is very good.

subtitle = sprintf("Forecasts made over %s to %s",
                   format(min(forecast_dates), "%B %d, %Y"),
                   format(max(forecast_dates), "%B %d, %Y"))

plot_canonical(results, x = "ahead", y = "ae", aggr = mean) +
  labs(subtitle = subtitle, xlab = "Weeks ahead", ylab = "Mean AE") +
  theme(legend.position = "bottom") + 
  scale_y_log10()

plot_canonical(results, x = "ahead", y = "wis", aggr = mean) +
  labs(subtitle = subtitle, xlab = "Weeks ahead", ylab = "Mean WIS") +
  theme(legend.position = "bottom") + 
  scale_y_log10()

plot_canonical(results, x = "ahead", y = "coverage_80", aggr = mean) +
  labs(subtitle = subtitle, xlab = "Weeks ahead", ylab = "Mean Coverage") +
  theme(legend.position = "bottom") + 
  coord_cartesian(ylim=c(0,1)) + geom_hline(yintercept = .8, color="black")

AE, WIS, and coverage by forecast date

Top line conclusions:

1. We crush it until December.
2. Quite a bit worse in Dec/January, though better than aardvark.
3. Coverage is generally good at all aheads.

theme_set(theme_bw())
plot_canonical(results, x = "forecast_date", y = "ae", aggr = mean,
               grp_vars = c("forecaster","ahead"), facet_cols = "ahead") +
  labs(subtitle = subtitle, xlab = "forecast date", ylab = "Mean AE") +
  theme(legend.position = "bottom") + 
  scale_y_log10()

plot_canonical(results, x = "forecast_date", y = "wis", aggr = mean,
               grp_vars = c("forecaster","ahead"), facet_cols = "ahead") +
  labs(subtitle = subtitle, xlab = "forecast date", ylab = "Mean WIS") +
  theme(legend.position = "bottom") + 
  scale_y_log10()

plot_canonical(results, x = "forecast_date", y = "coverage_80", aggr = mean,
               grp_vars = c("forecaster","ahead"), facet_cols = "ahead") +
  labs(subtitle = subtitle, xlab = "forecast date", ylab = "Mean Coverage") +
  theme(legend.position = "bottom") + 
  coord_cartesian(ylim=c(0,1)) + geom_hline(yintercept = .8, color="black")

Median relative WIS

Relative to baseline; scale first then take the median.

• Our November/December performance is the biggest issue. Though this is uncorrected data.

plot_canonical(results, x = "ahead", y = "wis", aggr = median,
               base_forecaster = "COVIDhub-baseline", scale_before_aggr = TRUE) +
  labs(subtitle = subtitle, xlab = "Weeks ahead", ylab = "Median relative WIS") +
  theme(legend.position = "bottom") + 
  geom_hline(yintercept = 1)

plot_canonical(results, x = "forecast_date", y = "wis", aggr = median,
               grp_vars = c("forecaster", "ahead"), facet_cols = "ahead",
               base_forecaster = "COVIDhub-baseline", scale_before_aggr = TRUE) +
  labs(subtitle = subtitle, xlab = "Forecast date", ylab = "Median relative WIS") +
  theme(legend.position = "bottom") + 
  geom_hline(yintercept = 1)

(Geometric) Mean relative WIS

Relative to baseline; scale first then take the geometric mean, ignoring a few 0’s. I think this is potentially more useful than the median/mean for relative WIS (or relative AE), but I haven’t completely thought it through. Putting the results here to be provocative.

geom_mean <- function(x) prod(x)^(1/length(x))
plot_canonical(results %>% filter(wis > 0), x = "ahead", y = "wis", 
               aggr = geom_mean,
               base_forecaster = "COVIDhub-baseline", scale_before_aggr = TRUE) + 
  labs(subtitle = subtitle, 
       xlab = "Weeks ahead", ylab = "Mean (geometric) relative WIS") +
  theme(legend.position = "bottom") + 
  geom_hline(yintercept = 1)

plot_canonical(results %>% filter(wis > 0), x = "forecast_date", y = "wis", 
               aggr = geom_mean, facet_cols = "ahead",
               grp_vars = c("forecaster", "ahead"),
               base_forecaster = "COVIDhub-baseline", scale_before_aggr = TRUE) +
  theme(legend.position = "bottom") + 
  labs(subtitle = subtitle, 
       xlab = "Forecast date", ylab = "Mean (geometric) relative WIS") +
  geom_hline(yintercept = 1)

Scores by target date (not forecast date)

• Again, our Achilles heel is November/December, but we’re using uncorrected data.

plot_canonical(results, x = "target_end_date", y = "wis", aggr = mean,
               dots = TRUE, grp_vars = "forecaster") + 
  labs(subtitle = subtitle, xlab = "Target date", ylab = "Mean WIS") +
  theme(legend.position = "bottom") + 
  scale_y_log10()

plot_canonical(results, x = "target_end_date", y = "wis", aggr = mean,
               dots = TRUE, grp_vars = c("forecaster", "ahead"), 
               facet_cols = "ahead") +
  labs(subtitle = subtitle, xlab = "Target date", ylab = "Mean WIS") +
  theme(legend.position = "bottom") + 
  scale_y_log10()

Maps (mean score over forecast dates and aheads)

maps <- results %>%
  group_by(geo_value, forecaster) %>%
  summarise(across(wis:ae, mean)) %>%
  pivot_longer(wis:ae, names_to = "score") %>%
  group_by(score) %>%
  mutate(time_value = Sys.Date(),
         r = max(value)) %>%
  group_by(forecaster, .add = TRUE) %>%
  group_split()


maps <- purrr::map(maps, ~as.covidcast_signal(
  .x, signal = .x$score[1], data_source = .x$forecaster[1], geo_type = "county"))
maps <- purrr::map(maps, 
                   ~plot(.x, choro_col = scales::viridis_pal()(3), 
                         range = c(0, .x$r[1])))
nfcasts <- length(unique(results$forecaster))

Mean AE

cowplot::plot_grid(plotlist = maps[1:nfcasts], ncol = 3)

Mean WIS

cowplot::plot_grid(plotlist = maps[(nfcasts+1):length(maps)], ncol = 3)

Trajectory plots

pd <- evalcast:::setup_plot_trajectory(
  bind_rows(qar_dc, submitted %>% filter(forecaster == "CMU-TimeSeries")), 
  intervals = 0.8,
  geo_type = "county", 
  start_day = min(forecast_dates) - 60)

g <- ggplot(pd$truth_df, mapping = aes(x = target_end_date))

# build the fan
g <- g + geom_ribbon(
  data = pd$quantiles_df,
  mapping = aes(ymin = lower, ymax = upper, fill = forecaster, 
                group = interaction(forecaster,forecast_date)),
  alpha = .1) +
  scale_fill_viridis_d(begin=.25, end=.75)

# line layer
g <- g +
  #geom_line(aes(y = .data$value.y), color = "#3182BD") + # corrected
  geom_line(aes(y = value)) + # reported
  geom_line(data = pd$points_df, 
            mapping = aes(y = value, color = forecaster, 
                          group = interaction(forecaster,forecast_date)),
            size = 1) +
  geom_point(aes(y = value)) + # reported gets dots
  geom_point(data = pd$points_df, 
             mapping = aes(y = value, color = forecaster),
             size = 3) +
  scale_color_viridis_d(begin=.25, end=.75)

g + theme_bw(base_size = 20) + 
  facet_wrap(~geo_value, scales = "free_y", ncol = 5) +
  theme(legend.position = "top") + ylab("") + xlab("")